###############################
# Analysis of conjoint data 
###############################

rm(list=ls())

# load packages 
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##############################
# Prepare data
##############################

# load data
d = read.dta("~/Dropbox/Cerrillos 2017/08_replication/conjoint_cerrillos_LAPS.dta")
dim(d)

# check data
d$idnum = as.numeric(d$idnum)
d$failcon[is.na(d$failcon)]=0
table(d$failcon, exclude = NULL) 
table(d$enumerator)
table(d$selected, exclude = NULL)

# drop failcon and non-responses
d = subset(d, d$selected!=99)
d = subset(d, d$selected!=88)
d = subset(d, d$failcon==0)
d = subset(d, d$failcon!="NA")
table(d$selected, exclude = NULL)
table(d$failcon, exclude = NULL)

##################
# Diagnostics
##################

# Balance check
prop.table(table(d$atideology))
prop.table(table(d$atprofession))
prop.table(table(d$atage))

# Profile order effect
d$fpair = as.factor(d$pair)
d1 <- lm(d$selected ~ d$atideology + d$atprofession + d$atage + 
                      d$fpair + 
                      d$atideology*d$fpair + d$atprofession*d$fpair + d$atage*d$fpair)
d1_c <- round(coeftest(d1, vcov = vcovCluster(d1, cluster = d$idnum)),4)
d1_c 

# Randomization female
table(d$b6)
d$female = 0
d$female[d$b6==1]=1
table(d$female)

d2 <- lm(d$female ~ d$atideology + d$atprofession + d$atage)
d2_c <- round(coeftest(d2, vcov = vcovCluster(d2, cluster = d$idnum)),4)
d2_c 

# Randomization enumerator
table(d$enumerator, exclude = NULL)

d3 <- lm(d$enumerator ~ d$atideology + d$atprofession + d$atage)
d3_c <- round(coeftest(d3, vcov = vcovCluster(d3, cluster = d$idnum)),4)
d3_c 

###################
# Appendix H
####################

d1_c
stargazer(d1_c,title="Profile order effects",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          keep = c(10:29),     
          covariate.labels=c("Left*Pair2",
                             "Left*Pair3",
                             "Left*Pair4",
                             "Left*Pair5",
                             "Professor*Pair2",
                             "Engineer*Pair2",
                             "Professor*Pair3",
                             "Engineer*Pair3",
                             "Professor*Pair4",
                             "Engineer*Pair4",
                             "Professor*Pair5",
                             "Engineer*Pair5",
                             "40*Pair2",
                             "50*Pair2",
                             "40*Pair3",
                             "50*Pair3",
                             "40*Pair4",
                             "50*Pair4",
                             "40*Pair5",
                             "50*Pair5"),     
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          #star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Dependent variable:",
          dep.var.labels   = "Electoral Choice",
          out="~/Dropbox/Cerrillos 2017/08_replication/t_profile_order.tex")     

d2_c
stargazer(d2_c,title="Balance check (female)",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          omit = c("Constant"),
          covariate.labels=c("Left",
                             "Teacher",
                             "Engineer",
                             "40",
                             "50"),     
          no.space=TRUE,  
          table.placement = "H",
          dep.var.caption  = "Dependent variable:",
          dep.var.labels   = "Female",
          out="~/Dropbox/Cerrillos 2017/08_replication/t_balance_female.tex")     

d3_c
stargazer(d3_c,title="Balance check (enumerator)",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          omit = c("Constant"),
          covariate.labels=c("Left",
                             "Teacher",
                             "Engineer",
                             "40",
                             "50"),     
          no.space=TRUE,  
          table.placement = "H",
          dep.var.caption  = "Dependent variable:",
          dep.var.labels   = "Female",
          out="~/Dropbox/Cerrillos 2017/08_replication/t_balance_enumerator.tex")     
